Relevance of visco-plastic theory in a multi-directional inhomoge- 
neous granular flow 



p. -P. CORTET^'2, D. BONAMY^, F. Daviaud^, O. Dauchot^, B. Dubrulle^ and M. Renouf^ 



o 
o 

> 
o 

:^ 



o 



7d 

o 
o 



1 CEA, IRAMIS, SPEC, CNRS URA 2464, Grp. Instabilites & Turbulence, 91191 Gif-sur-Yvette, France 

2 CEA, IRAMIS, SPCSI, Grp. Complex Systems & Fracture, 91191 Gif-sur-Yvette, France 

^ LaMCoS, CNRS UMR 5259, INS A Lyon, 18-20 rue des sciences, 69621 Villeurbanne, France 



PACS 45 . 70 . -n - Granular systems 

PACS 83 . 70 . Fn - Granular solids 

PACS 46 . 10 . +z - Mechanics of discrete systems 

Abstract. - We confront a recent visco-plastic description of dense granular flows [P. Jop et 
a/., Nature, 441 (2006) 727] with multi-directional inhomogeneous steady flows observed in non- 
smooth contact dynamics simulations of 2D half- filled rotating drums. Special attention is paid to 
check separately the two underlying fundamental statements into which the considered theory can 
be recast, namely (i) a single relation between the invariants of stress and strain rate tensors and 
(ii) the alignment between these tensors. Interestingly, the first prediction is fairly well verified over 
more than four decades of small strain rate, from the surface rapid flow to the quasi-static creep 
phase, where it is usually believed to fail because of jamming. On the other hand, the alignment 
between stress and strain rate tensors is shown to fail over the whole flow, what yields an apparent 
violation of the visco-plastic rheology when applied without care. In the quasi-static phase, the 
particularly large misalignment is conjectured to be related to transient dilatancy effects. 
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Despite granular flows ubiquity in nature and indus- 
try, their description still remains an open issue. Three 
regimes are usually distinguished: The rapid dilute flow 

- gaseous-like - regime, where grains interact through bi- 
nary collisions, is described within the framework of the ki- 
netic theory [1]; The slow flow - solid-like - regime, where 
only friction between grains is relevant, is described us- 
ing soil mechanics [2]. In between, there is a dense flow 

- liquid-like - regime, where inertia becomes important, 
friction remaining however relevant. This last regime has 
been widely investigated experimentally and numerically 
(see [3] for a review) in various flow configurations and has 
been the focus of many theoretical works [4-8]. However, 
while these descriptions reproduce some of the features ob- 
served experimentally, a unified description is still missing. 
In particular, the possibility for a single and local consti- 
tutive law within the framework of continuum mechanics 
is still debated. 

Progresses have been made recently in this context. Di- 
mensional analysis [3,9,10] was used to show that, in sim- 
ple incompressible unidirectional uniformly sheared flows, 
bead diameter grain density p, global shear stress r, 
global pressure P and global shear rate 7 should be re- 



lated as r/P = /i(/), where /i(/) is an effective fric- 
tion coefficient, function of the so-called Inertial Number 
I = "yd/ \fPj~p expected to control locally the rheology 
of the material. Such relation has first been evidenced, 
numerically, in plane shear [9, 10] and, experimentally, in 
inclined plane [3] configurations. Local tensorial extension 
of this relation was recently proposed by Jop et al. [11] as 
a constitutive law for dense granular flows: 
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2lijlij ^he norm of its deviatoric part 

itj = 7ij - l/27/c/c^ij, P = -l/2cFkk the local pressure 
and Vi the components of the velocity. In the compress- 
ible case, a consistent local rheology would imply in addi- 
tion a univocal dependance v = v{I) of the volume frac- 
tion on the inertial number [10]. Introducing equation ([T]) 
in the momentum equation, Jop et al. succeeded to fit 
the surface velocity profile for a steady unidirectional flow 
down an inclined plane with walls. More recently, rela- 



p-i 



P.-P. Cortet et al 



tion ([T]) was confronted with simulations of complex multi- 
directional non-steady inhomogeneous flow - the collapse 
of a cylinder of granular matter onto a horizontal plane [12] 
- and was shown to hold fairly well within the fast flow- 
ing region i.e. for / ranging from 10~^ to 0.8. In this 
last work, relation ([T]) has however been verified only for 
positions where stress and strain rate tensors had parallel 
principal directions, excluding 10% of data from the test. 
Therefore, despite these various experimental and numer- 
ical results, the relevance and limits of the visco-plastic 
relation ([T]) remains to be investigated in granular flows 
involving small values of inertial number / - in which the 
stress and strain rate tensors are particularly prone to be 
non-aligned - as e.g. avalanches flows where large scale 
non-local effects are expected [13-16]. 

In the present Letter, we confront this visco-plastic de- 
scription of granular flows with steady surface flows ob- 
served in a 2D rotating drum, using extensive numerical 
simulations providing large statistics. The motivation for 
studying this specific flow is double since, still remaining 
in a "simple" 2D geometry, (i) it provides a general highly 
multi-directional inhomogeneous situation and (ii) it offers 
a wide range of inertial number: 10~^ < / < 5 x 10~^. In 
such a flow, we report an apparent failure of the consti- 
tutive relation ([T]) over the whole range of /, when ap- 
plied in the "natural" frame of the flow. To overcome this 
failure, we highlight the fact that equation ([T]) is actu- 
ally composed of two independent results, namely (a) a 
single relation between the invariants of stress and strain 
rate tensors and (b) the alignment between these tensors, 
that constitute its two underlying fundamental hypothe- 
sis. The system of these two conditions actually recasts 
equation ([T]) in a more fundamental formulation. Then, 
the norms of the stress and strain rate tensors are found 
to be univocally related through the invariant relation de- 
rived from equation ([T]) over the whole range of /, meaning 
condition (a) is verified over the whole flow. Finally, the 
failure of the tensorial relation ([T]) is shown to result from a 
misalignment between the stress and strain rate tensors - 
meaning condition (b) fails - contrary to what is assumed 
in previous works. Possible origins of this misalignment 
are discussed: In the deep "solid" part of the flow, the total 
deformation experienced by the granular packing during 
its lifetime is less than one bead diameter. This induces 
transient dilatancy effects - from a Lagrangian point of 
view - that could be at the origin of the persistence of the 
observed large misalignment. 

The simulations have been performed using the Non- 
Smooth Contact Dynamics method [18]. The scheme has 
been described in detail in [19] and is briefly recalled 
below. An immobile drum of diameter Dq = 45cm is 
half-filled with 7183 rigid disks of density p = 2.7g.cm~^ 
and diameter uniformly distributed between 3 and 3.6mm. 
The weak polydispersity introduced in the packing pre- 
vents 2D ordering effects. The normal restitution coeffi- 
cient between two disks {resp. between disks and drum) 
is set to 0.46 {resp. 0.46) and the friction coefficient to 
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Fig. 1: (a-c) time and ensemble averaged streamwise velocity 
Vx for a Q = 6rpm rotation velocity: (a) map with normal- 
ized colormap, (b) with linear and (c) logarithmic z-axis as a 
function of the depth z at the center of the drum (a: = 0). (d) 
volume fraction zy as a function of the depth z at the center of 
the drum (x = 0). Velocities are non-dimensionalized by ^/gd 
where g is the gravity constant and d the mean bead diame- 
ter. In (b) and (c), solid lines are the result of a linear and an 
exponential data fit. 

0.4 {resp. 0.95). Once the initial packing is stabilized, 
a constant rotation velocity Q ranging from 2 to 12 rpm 
is imposed to the drum. After one round, a steady con- 
tinuous surface flow has developed. Then, one starts to 
capture 400 snapshots equally distributed over one rota- 
tion of the drum. For each bead, one records the position 
of its center of mass and its instantaneous velocity ~c 
over the time step dt = 6x 10~^s. For each rotation veloc- 
ity, we performed 20 experiments starting from different 
initial packings. The reference frame 5R is defined as the 
frame rotating at Q with the drum that coincides with the 
reference frame % = (e^, e^) (cf. fig. [l]), fixed in the 
laboratory, where {resp. ej) is parallel {resp. perpen- 
dicular) to the mean flow free surface. The mean values of 
the various Eulerian fields (volume fraction, velocity and 
stress tensor) are obtained, at the grain scale (see [19]) 
and without any further coarse-graining, after both time 
and ensemble averaging over all the 400 frames of all the 
20 experiments. 

Steady surface flows in rotating drums provide a con- 
tinuum of flow regimes from "quasi-static' to "liquid". 
Figures [ijb) and (c) presents typical velocity profiles 
Vx{x = 0, z) measured in our numerical simulations "at 
the center" of the drum {i.e. along the z- aligned slice at 
X = 0). Two phases can be identified: A flowing layer ex- 
hibiting a linear velocity profile and a "solid" phase expe- 
riencing creep motion where Vx decays exponentially with 
z. Such velocity profile corresponds well to experimental 
measurements in rotating drum [3,20-22] as well as in heap 
geometry [3,23]. It is worth mentioning that the charac- 
teristic length ^ associated with the exponential creep is 
found to be ^ = 12 ± 2 in our numerical drum over the 
range of explored rotating velocity. This is significantly 
larger than that reported by Crassous et al. [24] in heap 
geometry, ^ = 1 ± 0.2. This difference may emphasize the 
role played by the forcing on the selection of the veloc- 
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Fig. 2: Effective friction coefficient Txz\i^\hizP as a function 
of tfie inertial number for time and ensemble averaged data, 

(a) and (b) over tfie wfiole 2D drum, (c) and (d) at tfie center 
of tfie drum (x = 0), witfi linear, (a) and (c), and logarithmic, 

(b) and (d), x-axis. Green/light gray points correspond to 
Q = 12rpm, red/intermediate gray points to Q = 6rpm, and 
blue/dark gray points to O = 2rpm. 



ity profiles within the bulk of the material. Indeed, the 
forcing in rotating drum and heap are very different since, 
from a Lagrangian viewpoint, the gravity is rotating in 
the first case while it remains frozen in the second. The 
volume fraction is found to be almost constant, around 
0.827±0.05, in the creeping layer while it decreases slightly 
in the flowing layer, down to 0.80 at z = —4 (cf. fig. [ijd)) 
as expected from dilatancy effects. 

The full 2D fields of deviatoric stress r^j(x,z), devia- 
toric strain rate jfj{x^z) and pressure P{x^z) are then 
computed to test the relevance of constitutive equation 
([T]). Figures |2| a) and (b) show the variations of the ra- 
tio Txz\i^\/liz-P ^ function of / for time and ensemble 
averaged data corresponding to the whole 2D drum and 
for three different rotation velocities. If the rheology were 
correct, these data would collapse on the single master 
curve This is clearly not the case as we observe 

an important dispersion of the data - even for the larger 
values of /. Additionally, as one can see in figures [2jc) 
and (d), the collapse predicted by equation ([T]) even does 
not work well (though better) when restricting data to the 
slice X = located at the center of the drum, where the 
flow is expected to be locally unidirectional (/ to etc) and 
almost homogeneous {dx — 0) as in the plane shear con- 
figuration. Moreover, the shape of these curves does not 
correspond to the logarithmic shape previously evidenced 
in [3,10,11], especially, because of the unexpected bump 
of rxz\Y\hizP around / = 10"^ - lO""^. 

The multi-directionality of the considered flow implies 
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Fig. 3: Invariant effective friction coefficient \r\/P as a function 
of the inertial number for time and ensemble averaged data, (a) 
and (b) over the whole 2D drum, (c) and (d) at the center of 
the drum (x = 0), with linear, (a) and (c), and logarithmic, 
(b) and (d), x-axis. Same color/gray level code as in figure [2] 

that no reference frame appears more natural than an- 
other to study tensors. Therefore, it is reasonable to re- 
cast tensorial equations into reference frame invariant rela- 
tions. For granular flows, searching for a rheology consists 
mainly in relating the deviatoric stress tensor r^j to the 
deviatoric strain rate tensor jf,. Since both of these ten- 
sors are symmetric with zero trace, both of them are fully 
characterized through one invariant only, e.g. their norm, 
and their principal directions. Accordingly, the rheology 
of equation ([T]) can be recast in the following invariant 
way: 

(a) ^=mW, 

(b) Tij and jf, have identical principal directions. 



where |r| = y is the deviatoric stress tensor norm. 

In unidirectional flows, such as those observed in plane 
shear geometry [10] or down to rough inclines [11], con- 
dition (b) is naturally ensured by the flow symmetries so 
that verifying equation ([T]) reduces to verifying point (a) 
only. As we will see thereafter, this is no more the case 
for a multi-directional flow for which the alignment of the 
stress and strain rate tensors cannot be assumed [17,25]. 
The split of equation ([T]) into the two independent condi- 
tions (a) and (b) recasts the rheology of Jop et al. into its 
two underlying fundamental assumptions which are com- 
pletely independent. In the sequel, we will see that the in- 
variant form of the visco-plastic rheology - condition (a) 
- actually holds well over the whole rotating drum and 
the whole range of / whatever condition (b) holds or not. 
This reformulation of equation (IT]) actually allows to test 
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the visco-plastic theory for very small values of / down to 
10~^ where the condition on stress and strain rate tensor 
alignment is largely violated. 

Let us first consider point (a). In figure [sj we have plot- 
ted the invariant effective friction coefficient \t\/P as a 
function of the inertial number for the same data as in fig- 
ure|2] We see that all the data appear to collapse quite well 
on a single master curve independently of the rotation 
velocity Vt. The shape of this master curve and the value 
of the corresponding friction coefficient are perfectly con- 
sistent with the ones already evidenced (/i = 0.30±0.10 for 
/ < 0.1) in different granular unidirectional ffows [3,10,11]. 
It is important to emphasize that this collapse is a non 
trivial result since we have plotted data for the whole 2D 
highly inhomogeneous fields. In figures lajc) and (d), we 
restrict data to the slice x = located at the center of the 
drum where the fiow is expected to be almost "plane shear 
like". We see that the \t\/P profiles for the three consid- 
ered rotation velocities fall with a very good precision on 
the same master curve x = 0) for values of / as small 
as 10~^. Despite this last remark, this collapse is a strong 
evidence for the possible relevance of the constitutive law 
of equation ([T]) in steady inhomogeneous multi-directional 
fiows. More importantly, it is the first time such frictional 
rheology seems to hold over such a low inertial number 
range (10~^ < / < 5 x 10~^) whereas it is often expected 
to fail because of the jamming transition for / < 10~^. 
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Fig. 4: (a) Volume fraction v and (b) angle difference between 
principal directions of Tij and 7^^- as a function of the inertial 
number for time and ensemble averaged data corresponding to 
the whole 2D drum. Same color/gray level code as in figure |2] 

We already said that in compressive fiows a local rheol- 
ogy would imply a single dependence of the volume frac- 
tion V with the inertial number /. In figure [4]^a), we see 
that the data corresponding to the whole 2D fiow collapse 
indeed quite well on a master curve v = supporting 
partially the proposed rheological approach. 

The second step is now to check point (b) i.e. the 
collinearity between the deviatoric stress r^j and strain 
rate jfj tensors. Figure ^h) shows the angle difference 
56 between the principal directions of r^j and jfj as a 
function of the inertial number. The collinearity is not 
verified since 50 ranges from ±10°, for the larger iner- 
tial numbers (5 x 10~^ < / < 5 x 10~^) correspond- 
ing to the "liquid" surface fiow, up to 35° in the deeper 
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Fig. 5: (a) Map of the deviation \r\/P - \r\/P{x = 0,z) of 
the effective friction coefficient \r\/P from the one = 0) 

at the center of the drum for Q — 6rpm. (b) angle difference 
between principal directions of nj and jfj as a function of the 
inertial number for Q = 12rpm. In (b), data point colors are 
mapping the x position in the drum. 



part of the ffow, corresponding to smaller inertial numbers 
(10~^ < / < 5 X 10~^). Such increase in the dispersion 
of the tensors misalignment, as / decreases, has also been 
reported in shear zones [17]. More importantly, there is 
no univocal relation ruled by the inertial number such as 
50 = 50{I). Consequently, the fact the tensorial law ([T]) is 
not verified contrary to its invariant version, relation (a), 
on the norms of the tensors holds is the direct consequence 
of the misalignment of the stress and strain rate tensors. 
Moreover, it shows, looking at the larger inertial numbers 
in figures [2] and [4]^b), that a 10° discrepancy between prin- 
cipal directions is already significant to make the visco- 
plastic tensorial rheology, in its non-invariant version of 
equation ([T]), fiatter. 

In figure |3) the collapse of the data on a master curve 
is not perfect and one may wonder if the dispersion 
around the "mean" friction coefficient x = 0) mea- 
sured at the center of the drum is only due to a lack of 
statistics. Therefore, in figure [sja), we plot the map of 
the deviation of the effective friction coefficient from the 
one at the center of the drum. One can observe a spatial 
structuration with a quadripolar pattern. This spatial de- 
pendence of the profile \r\/P = jii{I^x) with the position 
X may be the consequence of non-local effects reffecting 
the long range inffuence of the boundaries through force- 
chains, but remains however an unexplained feature. It 
is worth to notice that a pure geometrical dependence of 
\r\/P has been reported in [25] where the interest of study- 
ing the present regime of mixed geometrical and shear rate 
dependence was highlighted. In parallel, in figures [5|b) 
and[6]^a), we see that the cloud of data 50 vs. I shows 
as well a non-trivial spatial structure. There is however 
no evidence for any link between the spatial structures of 
the angle 50 and of the deviation around = 0), the 

origins of which are still to be understood. 

Previous studies [17,25] have shown that the stress and 
the strain rate tensor principal directions can differ sig- 
nificantly due to transient dilatancy effects as long as the 
initial conditions of the packing have not been forgotten. 
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Fig. 6: (a) Map of the angle difference between principal direc- 
tions of Tij and jfj for Q — 12rpm. (b) Map of the logarithm of 
the norm of the strain rate tensor times the half drum rotation 
duration log \j^\/2Q for Q = 12rpm. We indicate with a black 
line the limit where this deformation 6 = |7^|/2f^ = 1. 



i.e. as long as the local deformation is less than ~ 1. 
Even if our flow is stationary, the transient effect reported 
above could be the explanation for the large misalignment 
observed in the deep part of the rotating drum (cf. fig. 
|6ja)). Indeed, in the "solid" part of the granular matter, 
a given "Lagrangian" packing is created at the bottom 
of the surface avalanche. Then, it has a short lifetime 
reduced to half a drum rotation period 1/2Q correspond- 
ing to the time needed for him to reach the top of the 
avalanche where it is destructed. During this process, the 
deformation experienced by the packing is quite small: in 
figure [6]^b), we plot the map of the logarithm of the norm 
of the strain rate tensor times the half drum rotation du- 
ration log 1 7^ 1/21] as well as a black line delineating where 
the deformation S = 17^^1/2^7 equals unity. The "solid" 
granular matter under this line is experiencing a transient 
deformation since it is less than one. Finally, we see that 
this "transient" zone matches quite well the zone (cf. fig. 
|6ja)) where \S0\ is larger than 10° making transient dila- 
tancy effects a relevant explanation for the persistence of 
large alignment discrepancies. This observation does not 
imply that the spatial structure observed for the angle 50 
is due to transient effects. On the contrary, this structure 
is created by the specific rotating drum boundary condi- 
tions, the small deformations experienced by the matter 
explaining then the fact that this structure of SO{x^ ^) 7^ 
can live on. 

Conclusion.-. — In this Letter, we have confronted 
a recent visco-plastic theory for dense granular flow to 
numerical steady flows in rotating drums. The existence 
of a univocal relation between a scalar local friction co- 
efficient - defined as the ratio between the norm of the 
deviatoric stress tensor to the pressure - and the so-called 
inertial number conjectured in the theory is in first ap- 
proximation fairly well verified over the whole drum, in 
both the creeping "solid" phase and the avalanching phase 
at the free surface, for inertial numbers ranging from 10~^ 
to 5 X 10~^. Small dilatancies evidenced at the flow sur- 
face are also found to be univocally related to the inertial 
number. On the other hand, the local alignment between 



deviatoric stress and strain rate tensors prescribed in the 
visco-plastic theory is found to fail significantly over the 
whole range of inertial number. This leads us to argue that 
the visco-plastic rheology proposed in [11] is extremely ef- 
ficient to describe quasi-unidirectional flows as those in- 
vestigated in [9-11] but remains unapplicable, when ap- 
plied in its non-invariant form (eq. ([T])), to highly multi- 
directional flows such as those observed in rotating drum. 
On the other hand, it is particularly remarkable that the 
Jop et al. rheology, when considered in its invariant form, 
holds for very small values of / down to 10~^ and that even 
if the stress and strain rate tensor misalignment is large. 
This result suggests in particular to reconsider the previ- 
ous apparent failure of the visco-plastic rheology observed 
for small / in other flow geometries [10,25,26] which could 
have been a partially wrong conclusion. The origin of the 
observed misalignment between tensors is conjectured to 
be the consequence of transient dilatancy effects linked to 
the small deformation experienced by the "solid" part of 
the granular matter when considered from a Lagrangian 
point of view. Finally, a quadripolar spatial structuration 
of the dispersion around the perfect collapse of the data 
on the visco-plastic invariant rheology is tentatively at- 
tributed to non-local effects set by boundary conditions. 
In this context, stress distribution in immobile granular 
packings was found to be well described by Linear Elas- 
ticity [27]. It will then be interesting to see whether the 
introduction of an additional elastic field describing the 
stress distribution in an equivalent "immobile" packing 
can explain the observed discrepancy. Work in this direc- 
tion is currently under progress. 

* * * 
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